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Abstract 

Calculations have been carried out to determine the neutron flux in- 
duced in the earth’s atmosphere by galactic protons and alpha particles at 
solar minimum for a geomagnetic latitude of 42°N. Neutron flux spectra in 
the energy range from ~ 10“ 8 to ** 10 5 MeV at various depths in the atmo- 
sphere were calculated using Monte Carlo and discrete ordinates methods, 
and various comparisons with experimental data are presented. The magnitude 
and shape of the calculated neutron- leakage spectrum at the particular lati- 
tude considered support the theory that the cosmic-ray-albedo-neutron-decay 
(CRAND) meachanism is the source of the protons and electrons trapped in 
the Van Allen belts. 
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INTRODUCTION 


Some of the neutrons produced in the earth ? s atmosphere by cosmic -ray 
bombardment escape and subsequently decay within the magnetosphere into pro- 
tons and electrons. It was initially proposed by Singer [1958] that this 
cosmic-ray-albedo-neutron-decay (GRAND) mechanism is the source of the protons 
and electrons trapped in the Van Allen radiation belts. Previous efforts (e.g., 
Dragt et al . [1966]; Hess and Killeen [1966]) to test this hypothesis have 
shown serious discrepancies between measured trapped-proton intensities and 
the intensities predicted by the CRAND theory. For example, Dragt et al . 

[1966] found the theoretical intensity at 55 MeV to be a factor of 50 too 
small. However, these predicted intensities have been based on neutron al- 
bedo spectra obtained from rather approximate calculations [Lingenfelter, 

1963] , especially for that portion of the albedo spectra £ 10 MeV, and at 
least part of the discrepancy may be ascribed to uncertainties in the mag- 
nitude of the neutron albedo spectra used in predicting the proton intensi- 
ties. Since the time of the Lingenfelter calculations, there has been con- 
siderable progress in the development of calculational methods for estimating 
the nucleon-meson cascade induced in matter by high-energy hadrons (e.g., 
Armstrong et al. [1972]), and therefore the present calculations were under- 
taken in an effort to better define the magnitude of the neutron albedo by 
using a calculational method that requires fewer approximations than previous- 
ly used. 

In the present work, the Monte Carlo code HETC [Chandler and Armstrong, 
1972] was used to compute the production and transport of protons, charged 
pions, and neutrons (> 12 MeV) due to protons and alpha particles incident 
on the top of the atmosphere. The neutron production £ 12 MeV from the HETC 
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calculation was used as input to the discrete ordinates code ANISN [Engle, 

1967] to obtain the low-energy (£ 12 MeV) portion of the neutron spectra. 

The basic input for this method of calculation is the proton and alpha-particle 
spectra incident at the top of the atmosphere. This differs from most methods 
used previously to calculate the neutron albedo (e.g., Lingenfelter [1963]; 
Newkirk [1963]; Hess et at, [1961]) in which the low-energy (<^ 10 MeV) neu- 
tron production was estimated from rather sparse experimental data, and only 
the low-energy 10 MeV) transport calculations were carried out in detail. 
Thus, the magnitude of the neutron fluxes computed here are based on experi- 
mental data for the galactic proton and alpha- particle spectra, whereas the 
normalization for the neutron fluxes from previous calculations has been 
determined from measurements of the low-energy neutrons in the atmosphere. 

In the next section the calculational details are given. The calcula- 
tions have been carried out for a geomagnetic latitude of 42°N and for solar- 
minimum proton and alpha-particle spectra. The calculations for incident 
alpha particles are more approximate than those for incident protons because 
in its present form the code HETC cannot treat alpha- particle -nucleus inter- 
actions, and therefore a rather crude model had to be used to obtain the 
description of the particles produced by alpha-particle interactions. 

In the last section the results are presented and discussed. Although 
the primary interest here is the neutron spectrum at the top of the atmosphere, 
the neutron spectra at various depths in the atmosphere have also been cal- 
culated and are compared with experimental data to provide a check on the 
accuracy of the calculational method for present purposes. 
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CALCULAT10NAL DETAILS 

Calculational Method: Incident Protons 

The Monte Carlo code HETC [Chandler and Armstrong, 1972] was used to 
compute the transport of protons, charged pions, primary alpha particles, 
and neutrons > 12 MeV, and the discrete ordinates code ANISN [Engle, 1967] 
was used for the low-energy (=£ 12 MeV) neutron transport. Since the standard 
version of HETC cannot transport alplia particles, the calculational method 
used for incident protons is discussed first, and the modifications to the 
code that were made to treat incident alpha particles in an approximate 
maimer are discussed later. 

Monte Carlo techniques, in conjunction with theoretical nuclear- inter- 
action models, are used in the high-energy code HETC to determine in detail 
the nucleon-meson cascade. The methods used in HETC for treating various 
physical processes are described in detail elsewhere [Armstrong et at . , 

1972], so only a brief description of those aspects of the code pertinent 
to the present problem will be described here. HETC takes into account 
charged-particle energy loss due to atomic ionization and excitation, ( non- 
elastic nucleon- nucleus and pion- nucleus collisions, elastic neutron- nucleus - 
collisions, and pion decay in flight and at rest. For nucleons above 3.5 GeV 
and pions above 2 . 5 GeV, nuclear interactions are treated using the intra- 
nuclear-cascade- extrapolation-evaporation model of Gabriel et al. [1970]. 

This model uses the particle-production data, obtained from an intranuclear- 
cascade calculation for intermediate- energy 3 GeV) nucleon- nucleus and 

pion- nucleus collisions, together with energy, angle, and multiplicity-scaling 
relations, which are consistent with the sparse e:xperimental data for high- 
energy collisions, to estimate the particle production for higher energy 
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(> 3 GeV) collisions. For nucleons between 12 MeV and 3.5 GeV and pions 
between 1.8 MeV and 2.5 GeV, nuclear interactions are treated using the 
intranuclear-cascade-evaporation model of Bertini and Guthrie [1971] • At 
each nuclear interaction that occurs during the Monte Carlo simulation of 
the nucleon-meson cascade, a calculation is performed to determine the energy, 
direction, and number of the interaction products, and the recoil energy, 
charge, and mass of the residual nucleus. The particles produced in nuclear 
interactions may be protons, neutrons, charged pions, and neutral pions from 
the intranuclear cascade and protons, neutrons, deuterons, tritons, 3 He 1 s , 
and alpha particles from the evaporation. The produced neutral pions are 
assumed to decay immediately into two photons, and these photons, as well as 
all photons produced in nuclear interactions, are neglected. The heavy par- 
ticles (A > 1) produced by evaporation are assumed to slow down and come to 

rest at their point of origin without undergoing nuclear interaction. Neu- 

* 

trons and protons above 12 MeV and charged pions above 1.8 MeV are followed 
until they eventually escape from the atmosphere, undergo nuclear absorption, 
or, in the case of charged pions, decay. Protons that slow down to 12 MeV 
or are produced below 12 MeV and charged pions that slow down to 1.8 MeV or 
are produced below 1.8 MeV are assumed to come to rest without undergoing 
nuclear interaction. Because of the relatively low density of the atmo- 
sphere, it is much more probable that negatively charged pions that come to 
rest will decay rather than undergo nuclear capture, so both positive and 
negatively charged pions reaching 1.8 MeV were assumed to decay. Muons from 
pion decay were neglected. The treatment of neutrons produced below 12 MeV 
is discussed later . 
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The accuracy of the code HETC for predicting the characteristics of 
nucleon-meson cascades has been checked for source -particle energies up to 
~ 30 GeV, and, in general, good agreement with experimental data has been 
found [Armstrong et al . , 1972]. Also, there is considerable indirect evi- 
dence that the code is reasonably accurate at energies much higher than 
30 GeV (e.g., Armstrong et al. [1972]; Armstrong and Alsmiller [1971]). 

Since the nuclear- interaction models used in HETC are inapplicable at 
very low energies, the neutron production £ 12 MeV computed by HETC was used 
as input to the discrete -ordinates code ANISN [Engle, 1967] to obtain the 
low-energy neutron transport. The neutron cross sections for the ANISN cal- 
culations were obtained from the ENDF/B cross section library [ENDF/B, 1970] 
and put in multigroup form by using the code XLACS [Greene et al . , 1972] . 

The ANISN calculations were made using 57 spatial intervals, S 16 P 3 angular 
quadrature, and 32 energy groups, with upscattering allowed in the lowest 
13 energy groups. Upscattering allows neutrons to gain as well as lose 
energy in elastic collisions with nuclei, and it is needed to properly pre- 
dict the shape of the neutron spectrum at very low (near thermal) energies . 
Calculational Method: Incident Alpha Particles 

It was necessary to resort to a very approximate method to treat inci- 
dent alpha particles. The position in the atmosphere at which an alpha par- 
ticle experienced its first nuclear collision was determined using a mean 
free path based on the cross section a. = 7r(r + r-) 2 , where r and r- are 
effective radii for fhe alpha particle and the struck nucleus (nitrogen and 
oxygen), respectively, and are given by r = 1.17 * 10“ 13 A 1 ^ 3 cm [Webber, 

1967]. The energy of the alpha particle at the collision site was obtained 

/ 

using the stopping-power formula for protons, with an approximate density- 
effect correction [Armstrong and Alsmiller, 1970] and scaling relations 
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(e.g., Barkas and Berger [1964]). To obtain the energy and direction of 
the collision products, a very approximate model [Gabriel, Santoro, and 
Alsmiller, 1971] was used in which it is assumed that the nucleons of the 
alpha particle enter the nucleus separately and independently except for their 
relative spatial locations when they enter the nucleus with each nucleon hav- 
ing a kinetic energy equal to one-quarter of the difference between the ki- 
netic energy and the binding energy of the alpha particle. From the products 
produced by all four of the incident nucleons and from the energy and momen- 
tum conservation, the charge, mass, and excitation energy of the residual 
nucleus were determined, and an evaporation calculation for the collision was 
then made. These assumptions allow the products from alpha-particle colli- 
sions to be computed using the standard collision models contained in HETC. 

The cascade initiated by the alpha-particle collision products was calculated 
in exactly the same manner as described previously for incident protons . 

Other Calculational Details 

The calculations have been made for an isotropic flux of galactic pro- 
tons and alpha particles incident on an infinite slab of atmosphere (79 atom 
% N, 21% 0) 1033 g/cm 2 in thickness, and the ground was represented by a 
layer of Si0 2 ~ 500 g/cm 2 in thickness with a density of 1.8 g/cm 3 . The 
density variation of the atmosphere [U.S. Standard Atmosphere, 1962] was 
taken into account since the pion decay probability per unit distance is 
density- dependent . 

The temperature variation with altitude was taken into account in the 
generation of the low-energy (< 12 MeV) neutron cross sections for the ANISN 
code. The following values [Hess, Canfield, and Lingenfelter, 1961] for the 
temperature in various depth intervals (measured from the top of the atmo- 
sphere) were used: 
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0 - 

235 g/cm 2 : 

219°K 

235 - 

485 g/cm 2 : 

234°K 

485 - 

710 g/cm 2 : 

260°K 

710 - 

1033 g/cm 2 : 

280°K 

> 

1033 g/cm 2 : 

300°K 


The proton and alpha-particle spectra presented in the review article 
by McDonald [1969] were used. These spectra are based on several balloon 
and satellite measurements made near the solar minimum in 1965 with a power 
law in energy extrapolation at high O 10 GeV/nucleon) energies. 

The calculations have been made for a vertical cutoff rigidity of 
4.6 GV, which corresponds to a vertical cutoff energy of 3.8 GeV for pro- 
tons and 6.3 GeV for alpha particles, and (for a dipole afield) a geomag- 
netic latitude of \ = 42°. Incident protons in the energy range from 3.8 GeV 
to 200 GeV and incident alpha particles in the energy range from 6,3 GeV to 
800 GeV were considered in the calculations. The geomagnetic cutoff at non- 
vertical angles was taken into account although at the latitude considered 
here the cutoff is rather sharp; i.e. , no particles with energies less than 

" f • . " / ' 

the vertical cutoff energy penetrate the magnetic field and essentially all 
particles with energies greater than the vertical cutoff are transmitted. 

The details of the procedure used in selecting source-particle energies > and 
directions with the effect of the geomagnetic cutoff included are discussed 

in Appendix A. . 
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RESULTS 

Neutron Flux Spectra vs Depth 

The calculated omnidirectional (over 4 tt) neutron- flux spectra at var- 
ious depths in the atmosphere are shown in Figs. 1 and 2. At depths of 
0 g/cm 2 (top of atmosphere), 200 g/cm 2 , and 1033 g/cm 2 (air-ground inter- 
face), calculated neutron spectra due only to incident protons are also 
shown to indicate the contribution of incident alpha particles to the neu- 
tron spectra at these depths. Error bars are shown for those, portions of 
the spectra (> 12 MeV) calculated using the Monte Carlo methods. The error 
bars shown in these and subsequent figures correspond to estimated statis- 
tical errors of one standard deviation. Statistical fluctuations are also 
associated with the results < 12 MeV since the neutron source distribution 
for the discrete ordinates calculation is obtained by the Monte Carlo methods. 
The statistical error < 12 MeV is estimated to lie generally between 5 and 
25%. 

The neutron spectrum calculated by Lingenfelter [1963] at 0 g/cm 2 is 
also shown in Fig. 1. This spectrum is referred to as the neutron flux 
spectrum in the Lingenfelter paper, but, according to commonly accepted 
nomenclature (e.g., Beckurts and Wirtz [1964]), the quantity actually cal- 
culated was the neutron current spectrum [Lingenfelter, 1972] . The neutron 
current spectrum at 0 g/cm 2 is also available from the present calculations 
and is shown in Fig. 1. The neutron spectrum calculated here is in good 
agreement with that calculated by Lingenfelter at low (< 10 MeV) energies, 
but the spectrum calculated here is considerably higher for energies ** 10 MeV. 
However, from the nature of the calculational method used by Lingenfelter, 
his spectrum is expected to be very approximate at energies £ 10 MeV. 


t 
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sphere. 
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1. Neutron spectra at various depths from the top of the atmo- 
Solar minimum, X ~ 42°N, - 
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The neutron- flux spectra measured in 1956-1957 by Hess et al. [1959] 
at depths of 200 and 1033 g/cm 2 are also shown in Fig. 2. The measured and 
calculated spectra differ somewhat at low energies, especially in the thermal - 
energy region, but are in quite good agreement at high energies. The ac- 
curacy of the measured spectra at very low energies («* 4 x 10’ 8 MeV) is 
estimated by Hess et at . [1959] to be about a factor of two. 

Neutron Flux vs Dept h 

I 

The depth dependence of the neutron flux in several energy ranges is 
shown in Fig. 3. The calculated fluxes are compared with the measurements 
of Boella et al . [1965] (A = 42°N, 1963) in the energy range from thermal 

to 20 MeV and with the measurements of Holt et al. [1966] (A = 42°N, 1964) 
in the 1- to 10-MeV region. The statistical error (one standard deviation) 

for the calculated results is generally less than 20%, with the greatest 
uncertainty at the larger depths. The fluxes show considerable variation 
near the air- ground interface. This variation is, of course, dependent upon 
the assumed composition for the ground. In the present calculations, Si0 2 
with a density of 1.8 g/cm 3 was used to simulate dry soil. The addition of 
water to the soil would have a significant effect on the low-energy neutron 
flux (see, for example, Yamashita et al. [1966]) because of the large ab- 
sorption cross section of hydrogen at low energies. Yamashita et al. [1966] 
measured the neutron flux over dry soil at sea level and obtained a flux 
about a factor of three smaller than that calculated here. Yamashita et al. 
[1966] point out that their measured neutron production at sea level is 
also about a factor of three smaller than the neutron production predicted 
from the data of Hess et al. [1959, 1961) and Newkirk [1963]. 0*Brien 
[1971] calculated the depth dependence of the neutron flux in the atmosphere 
using an essentially analytical propagation theory and a phenomenological 
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Fig. 2. Neutron spectra at various depths from the top of the atmo- 
sphere. Solar minimum, X = 42°N. 
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model for hadron-nucleus collisions. His calculations, which neglect the 
ground, predict neutron fluxes near sea level substantially lower than the 
fluxes calculated here. 

A comparison of the spatial dependence of the neutron flux in the 1- to 
10-MeV range near the top of the atmosphere at A » 42°N is shown in Fig. 4. 

The fluxes from the present calculations with and without incident alpha par- 
ticles included are shown, which indicate that for this particular latitude 
and energy and spatial regions incident alpha particles contribute « 171. The 
range of values shown in Fig. 4 from Haymes [1964] represent approximately the 
spread in the data obtained from four balloon flights. These measurements were 
made in 1963 with a detector sensitive to neutrons in the 1- to 14-MeV range. 
(According to the energy dependence of the flux measured by Haymes [1964] , the 
neutron flux in the 1- to 14-MeV range is approximately 10% higher than the 1- 
to 10-MeV flux* ) The results of the diffusion- theory calculations of Lingen- 
felter [1963] and the discrete ordinates calculations of Newkirk [1963] shown 
in Fig. 4 were taken from the paper of Haymes [1964]. (Lingenfelter* s results 
in this case are in terms of flux, not current [Lingenfelter, 1972].) Newkirk’s 
calculated results for A « 57°N and (presumably) 1961 were converted to A ~ 42°N 
by Haymes [1964] using the latitude dependence calculated by Lingenfelter [1963] . 
The measurements of Holt et at. [1966] were made at A = 42°N during the latter 
part of 1964 using a detector sensitive to neutrons in the 1- to 10-MeV range. 
Thus, the results shown in Fig. 4 are for cornnon conditions with the exceptions 
that the Haymes data are for 1 to 14 MeV instead of 1 to 10 MeV and the Haymes 
and Newkirk data are not at solar minimum. According to the measurements of 
Marker et at . [1969] , the 1- to 10-MeV neutron flux at A - 42°N can vary by 
« 201 during a solar cycle. 
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Fig. 4. Depth dependence of the flux of neutrons with energies between 
1 and 10 MeV near the top of the atmosphere. 
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Haymes [1964] and Holt et al. [1966] extrapolated their data near the 
top of the atmosphere to obtain the flux at 0 g/cm 2 , and these extrapolated 
fluxes are given in Table 1, together with calculated values for the flux 
at 0 g/cm 2 . Thus, at X » 42°N both the calculations and measurements give 
quite consistent results for the neutron flux in the 1- to 10-MeV region at 
the top of the atmosphere. 

Neutron Spectrum Near the Top of the Atmosphere 

In Fig. 5 portions of the calculated flux spectra at depths of 0 and 
5 g/cm 2 are shown and compared with the calculations of Lingenfelter [1963] , 
the measurements of Holt et al . [1966] and Preszler et al* [1972], and an 
approximate estimate made by Freden and White [1962]. At 5 g/cm 2 both the 
calculated and measured [Preszler et al . , 1972] spectra are for only the up- 
ward moving neutrons, i.e., the flux per unit energy due to neutrons having 
directions between 0 = 0° and 0 = 90° with respect to the zenith direction. 

As pointed out earlier, the Lingenfelter energy spectrum (solar mini- 
mum, X * 40°N) is for the neutron current, not the neutron flux. The Holt 
et al. [1966] spectrum is from data taken in late 1964 at X = 42 U N and has 
an E" 1,0 energy dependence. The measurements of Preszler et al. [1972] 
were made in 1971 at X ■ 40 o N. Freden and White [1962] estimated the neu- 
tron albedo spectrum needed to produce trapped proton spectra in agreement 
with proton measurements made in the radiation belt at L = 1.30 and B = 0.200 

o • 

Their prediction used the GRAND theory of injection and took into account 
proton losses in the atmosphere due to ionization and nuclear collisions 
but did not take into account the injection coefficient. Preszler et at. 
[1972] multiplied the Freden and White [1962] spectrum by seven to take into 
account the injection coefficient calculated by Dragt, Austin, and White 


TABLE 1 

1- to 10-MeV Neutron Flux at 0 g/cm 2 
(X « 42°N) 


Source 

m — — ww w i n wi 

Present calculations 
Calculated, Lingenfelter a 
Calculated, Newkirk 3 
Measured, Holt et at, [1966] 
Measured, Haymes [1964] 
a. From the Haymes [1964] paper. 


Flux 

(neutrons/cm 2 «sec) 
0.23 + 0.02 
0*26 
0,20 
0,2(3 

0,24 ± 0,02 
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Fig. 5. Neutron spectra near the top of the atmosphere. At S g/an 2 , 
both the calculated and measured spectra are for inward-moving neutrons. 


20 


[1966] , and it is the Freden and White spectrum as modified by Preszler et 
at. that is shown in Fig. 5. Leavitt [1972] recently measured the neutron 
albedo spectrum in the energy range from « 50 to « 150 MeV, and, although 
final analysis of the data has not yet been completed, Leavitt’s preliminary 
data indicate a spectrum considerably lower in magnitude than that measured 
by Preszler et at. [1972] . Zobel et at. [1972] measured the neutron-flux 
spectrum between about 3 and 60 MeV at 9 g/ cm 2 . Their measured spectrum 
over most of the energy range is lower by a factor of two or more than the 
spectrum at 10 g/cm 2 calculated here (Fig. 1). 

A comparison at a depth of 5 g/cm 2 of the omnidirectional flux spectrum 
and the calculated and measured [Preszler et at. 1972] flux spectra due to 
upward moving neutrons is shown in Fig. 6. White et at. [1972] reanalyzed 
the data of Preszler et at . to obtain the upward moving current spectrum, 
and Fig. 6 shows the measured and calculated current spectra due to upward 
moving neutrons and the calculated omnidirectional current spectrum. The 
spectra shown in Fig. 6 are defined as: 

4 > 4 7r (E) = the omnidirectional flux spectrum 

“ the angular flux spectrum integrated over 4 tt solid angle 

88 / 2TT dip / dy <f>(E,i p,v) 

0 -1 

^2tt® " the upward moving flux spectrum 

a the angular flux spectrum integrated over 2n solid angle 
about the zenith direction 

,2ir ' 1 , ' . ■■ v ‘ . V '/ 

■ / dip f du <P(E,\p,v) ; 

o o 


NEUTRON FLUX OR CURRENT PER UNIT ENERGY [neutrons /(cm 2 * sec ♦ MeV)] 
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Fig. 6. Neutron flux and current spectra at 5 g/cm 2 . The measured 
upward moving flux and current correspond to ^(E) and ^(E), respectively 



J^CE) = the omnidirectional current spectrum 

,27r , ,+i 

/ dip J dy | y | <|>(E,iJ>,y) ; 

0 -1 

J^Tf (E) » the upward moving current spectrum 

,2t r .1 

= / dip J dy y ij>(E,ip,y) ; 

0 0 

where 

<t>(E,ip,v) a the neutron flux per unit energy and solid angle at 
5 g/cm 2 depth, 
y = cose, 

e = the angle between neutron direction and zenith angle, 
ip = the azimuthal angle. 

According to the notation used here, Preszler et at, [1972] measured ^(E) 
and White et al. [1972] measured J2 7r (E). Figure 6 shows that in the energy 
region from 10 to 100 MeV about 801 of the omnidirectional flux (and cur- 
rent) at 5 g/cm 2 is due to upward moving neutrons. 

The contribution of incident alpha particles to the omnidirectional 
flux at 5 g/cm 2 is shown in, Fig. 7. At this depth, incident alpha particles 
contribute approximately 20, 35, and 45% to the flux spectrum in the energy 
ranges of 1 to 10 MeV, 10 to 100 MeV, and > 100 MeV, respectively. 






PRESENT CALCULATIONS (5 g/cnrT) -22 

PROTONS AND ALPHA PARTICLES — 

INCIDENT ~ 

PROTONS INCIDENT ~ 

MEASURED (5 g/cm 2 ), PRESZLER eta/. (1972) 
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DISCUSSION 

Dragt, Austin, and White [1966] used the neutron albedo spectrum cal- 
culated by Lingenfelter [1963] (shown in Fig. 5) and the theory of cosmic- 
ray- albedo -neutron- decay (CRAND) injection to arrive at a theoretical esti- 
mate of the proton intensity in the earth f s radiation belt. The theoreti- 
cal value of the proton flux at 55 MeV was found to be a factor of 50 smaller 
than the proton flux measured by Filz and Hoieman [1965] . Major uncertain- 
ties in the theoretical estimate were the albedo neutron flux and the mean 
atmospheric densities encountered by trapped protons, and Dragt, Austin, 
and White [1966] concluded that agreement between theory and experiment could 
be obtained only if the ratio of the albedo -neutron flux to the atmospheric 
density used in the calculations was increased by a factor of 50. Figure 5 
shows that the neutron leakage flux at 55 MeV calculated here is about a 
factor of 33 higher than the neutron leakage predicted by Lingenfelter [1963] . 
At a depth of 5 g/cm 2 , the measurements of Preszler et al. [1972] and the 
present calculations give fluxes at 55 MeV that are factors of about 45 and 
35, respectively, higher than the Lingenfelter [1963] curve. Thus, both the 
present calculations and the recent experimental data of Preszler et al. 

[1972] indicate that the CRAND theory is capable of correctly predicting the 
high-energy (£ 20 MeV) proton intensities observed in the inner zone of the 

Van Allen radiation belt. Neutron leakage spectra at other geomagnetic 

« 

latitudes are needed to allow revised injection [Dragt, Austin, and White, 
1966] and diffusion [Farley, Tomassian, and Walt, 1970] calculations to be 
made to obtain additional tests of the CRAND theory. 
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APPENDIX A 

SELECTION OF SOURCE -PARTICLE ENERGIES AND DIRECTIONS 
AND THE INFLUENCE OF THE GEOMAGNETIC FIELD 

First we consider the straightforward case of selecting energies and 
directions from the galactic proton or alpha-particle spectrum without re- 
gard to geomagnetic- field effects or any importance sampling. The omni- 
directional (over 4 tt) flux spectrum outside the geomagnetic field, <j> o (E) , 
is considered known and assumed to be isotropic; i.e., 


- 

<|>(E,fi) - 9 

where S is a unit vector along the direction of travel of the particle, 
corresponding current spectrum incident toward the earth is 


The 


J(E,ft) = <j>(E,fi) S*n , when 2s 0 , 

= 0 , otherwise, 

- .I. — _ - nAinnal f a arirl intn tllP 

WilCIC Jl JS.3 tile UUJLt V ev tVJ. UV/JUI 1 UJ. uv!, wtu U4.i Cv vv«* 


s»-hmnc:nVie' r ’e < 


The current spectrum can be written in terms of the properly nomalized 
source distribution, S(E,fi), for selecting energies and directions as 


where 


J(E,fi) = S(E,S) , 


$ f^rnax 
1 


dE / dfi J(E,5) 


E . 2tt 
min 




®min 


S(E,Q) = J(E,n)/[/ 


max / dn J(E,2)] , 

2tt 


E . 
nun 
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and E min and define the source-particle energy range being considered. 
By rewriting S(E,S), 


$ 

r*_(B) 1 

0 - dE 

Idty 

? 

<2> 

|2¥ 


[2p dy] 


- W Q Pl (E) dE p 2 W d* p 3 (y) dy , 


where 

ijj = the azimuthal angle, 
e = cos" 1 y * cos’ifnni) = zenith angle, 
W Q - 0/4, 

Pl (E) - 4 0 (E)/*, 

P 2 ft) = 1/2^, 
p 3 (y) - 2y. 


Source-particle energies and directions can then be realized by selecting 
E, t|>, and y according to the probability density functions (p.d.f.’s) Pj , 
p 2 , and p 3 , and giving each source particle the statistical weight 

Because the flux spectrum 4> 0 (E) decreases very rapidly with increasing 
energy, it is almost imperative that some importance sampling be used in 
the energy selection to obtain adequate statistics for the effects produced 
by high-energy source particles. Therefore, the following biased p.d.f. for 
E was used instead of p x (E) : 

p*(E) dE = p! p"(E|i) dE 


= Pi 


dE 


f 


where 
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p| « the fraction of source particles to be selected in energy 
group i, 

p**(E|i) * the p.d.f. for selecting E given that E lies within energy 
group i, 


and 



p|(E) dE 8 1 • 


Values for were assigned to give greater emphasis to the high-energy 
portion of the spectrum than specified by p x (E). The statistical weight 
for a source particle selected within energy group i then becomes 

W 8 W 0 Pl (E)/p|(E) 

- CE i+1 - E^/p! . 

Using the above procedure, the expected number of source particles in each 
batch selected in energy group i is N, where N is the total number 

of source particles per batch. The variance associated with ( n i ) was elim- 
inated by using quota sampling; i.e. , at the beginning of each batch <n^> 
source particles were assigned a priori to belong in the ith energy group. 

(M and p| were fixed so that <n. > would be an integral number.) 

The influence of the earth 1 s magnetic field can be taken into account 
by simply selecting E, 8, and ip in the manner described above and accepting 
or rejecting the selected coordinates according to the value of the cutoff 
energy at zenith angle 8 and latitude X. The cutoff energy is given by 
(e.g. , Hopper [1964]; Haffher [1967]) 


/ 

fZ P cX (9)1 

2 

/ 1 + 


- 1 

— 

me 2 i 



E cA ( 8) 8 me 2 
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where 

P c ^(0) = the cutoff rigidity 

4P v > 

r 2 (l + A - cos 3 A cose) 2 

P - the vertical cutoff rigidity 

V 

=14.9 cos 4 A (for dipole field), in units of GV , 

Z,mc 2 = the charge number and rest energy of the particle, 
r = (Rg + h)/Rg , 

Rg = the radius of the earth, 
h = the altitude. 

In the present calculations, an r of unity was used. 

It is of interest to note that the flux spectrum outside the geomag- 
netic field and the flux spectrum transmitted through the geomagnetic field 
can be simply related as 

♦ fl CE) 

♦JCH) = -^y- F x CH) , 

where 

= the omnidirectional (over 2n) flux at energy E and latitude A 
transmitted through the geomagnetic field and incident on 
the earth’s atmosphere, 

<j> o (E) = the omnidirectional (over 4w) flux at energy E outside the 
geomagnetic field, 

(E) « the flux transmission factor at energy E and latitude A. 

The expression for (E) can be arrived at by averaging the cutoff energy 
over all zenith angles; i.e.. 
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$f(E) = / HE, a) h.(Q) d 5 
A 4ir A 

<f>(E) »/* 

= — ^ — / h x (e) SAne 9 

o 

where 


h x (6) - 1 if B ^E cX (e) , 
» 0, otherwise. 


Thus, 


F X (E) 


/ c sine de 
0 


= 1 - A - sin 2 e c , 


where 

e c » the angle at which E ® E c , 

sin 2 e = -^Sl , 
c cos^x 

and 

4 P Z 
? _ v 

sasa ===== = = • 

M E/mc 2 + 1) 

Figures A1 and A2 show $*(E) at various X for protons and alpha particles, 
respectively, and support the statement made in the text that at the lati- 
tude of interest here (42°) the transmitted flux spectrum has a sharp energy 
cutoff. 

It is also of interest to consider the latitude dependence of the number 
of protons and alpha particles transmitted. Figure A3 shows the integral 
proton flux, T> , the integral alpha-particle flux, , and the ratio 

* A,p A, Cl 

r / * n , 

a/p A, a X,p 


30 









33 



Fig. A3. Integral proton flux, integral alpha-particle flux, and the 
alpha-particle- to -proton- integral -f lux ratio vs geomagnetic latitude at 
solar minimum. 
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